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Abstract 



The finite element computation of structures such as waveguides can lead to heavy 
computations when the length of the structure is large compared to the wavelength. 
Such waveguides can in fact be seen as one-dimensional periodic structures. In 
this paper a simple recursive method is presented to compute the global dynamic 
stiffness matrix of finite periodic structures. This allows to get frequency response 
functions with a small amount of computations. Examples are presented to show 
that the computing time is of order log2 N where is the number of periods of the 
waveguide. 

Key words: periodic structure, waveguide, finite element, recursive, dynamic, 
vibration 
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1 Introduction 



We study here the computation of structures considered as waveguides, as 
shown in figure 1, with symmetries which can be a translation (a), a rotation 
(b) or a periodicity (c). Thus these waveguides can be uniform or periodic. 
The vibration of such waveguides has been the topic of much research. One 
can find analytical or finite element models of waveguides and people are 
generally interested by the computation of wave propagations and dispersion 
curves or by the determination of the frequency response functions. A first 
approach considers structures with constant cross sections as the cases a) 
and b) of figure 1. For example, [1] and [2] used a wave approach to study 
the vibrations of structural networks composed of simple uniform beams, and 
solved for the dynamics of individual elements and of the junctions between 
elements by analj^ical methods. The efficiency is greatly improved compared 
to FE methods as a beam can be modelled using only a single element. 

The first numerical approaches were proposed by [3,4] to approximate the 
cross-sectional deformations by finite elements. The authors of [5,6] applied 
similar ideas to the calculation of wave propagations in rails using a finite 
element model of the cross-section of a rail. They then calculated dispersion 
relations and accelerances. Dispersion relations for elastic waves in helical 
waveguides were also considered by [7]. For general waveguides with a com- 
plex cross-section, the displacements in the cross-section can be described by 
the finite element method while the variation along the axis of symmetry is ex- 
pressed as a wave function. Following these ideas, [8,9,10,11,12,13] developed 
the spectral finite element approach. This leads to efficient computations of 
dispersion relations and transfer functions but special elements need to be de- 
veloped for each element type. This makes the connection with the standard 
use of the finite element method difficult and does not allow the benefits of 
powerful existing finite element software to be exploited. Similar techniques 
were also developed by [14,15] for the computation of dispersion relations in 
damped waveguides. 

More general waveguides can be studied by considering periodic structures. 
Numerous works provided interesting theoretical insights in the behaviour of 

these structures, see for instance the work of [16] and the review paper by 
[17]. Mead also presented a general theory for wave propagation in periodic 
systems in [18,19,20]. He showed that the solution can be decomposed into an 
equal number of positive and negative-going waves. The approach is mainly 
based on Floquet's principle or the transfer matrix and the objective it to 
compute propagation constants relating the forces and displacements on the 
two sides of a cell (a single period) and the waves associated to these con- 
stants. For complex structures FE models are used for the computation of the 
propagation constants and waves. The final objective is to compute dispersion 
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relations to use them in energetic methods, see [21,22,23,24,25]. In [26] the 
general dynamic stiffness matrix for a periodic structure was found from the 
propagation constants and waves. It leads to a matrix linking the extreme sides 
of the structure and allows to compute transfer functions in the structure. 

The last approach is purely computational and uses rotational and cyclic sym- 
metries to solve the problem by a decomposition of the displacements in cosine 
and sine functions. It is thus possible to find the transfer functions or modal 
shapes for periodic structures. A review of the current practises can be found 
in [27]. These methods allow computing the frequency response functions in 
a number of operations proportional to the number of cells in the structure. 
This paper develops this last approach and presents a recursive method to 
calculate the forced response of structures such as those illustrated in figure 1. 
A section of the waveguide is modelled using conventional FE methods, using 
a commercial FE package. The resulting mass, stiffness and damping matrices 
are then post-processed to give the dynamic stiffness matrix of the cell. Then 
a recursive method is applied to compute the global dynamic stiffness matrix 
of the whole waveguide and finally the transfer functions in the structure. 

This paper presents a different approach from the previously published pa- 
per [26]. Both papers aim at computing the global dynamic stiffness matrix 
of a N cells structure. But in [26] waves in a period were computed and from 
these waves the dynamic stiffness matrix of a complete structure was obtained 
with a computational cost independent of the number of periods in the struc- 
ture. However the computation of the waves can be time consuming when the 
number of dofs in a section is large because this needs the computation of 
eigenvalues of non symmetric matrices. On the contrary, in the present ap- 
proach no wave needs to be computed and the global dynamic stiffness matrix 
is obtained by products and inverses of matrices with the same dimensions as 
the dynamic stiffness matrix of a cell. Then, the frequency response functions 
can be obtained easily without the computation of any wave. 

The paper is divided into two parts. In the first part the recursive method for 
the finite element analysis of periodic structures is presented. In the second 
part two examples consisting in a beam and a plate are described before the 
conclusion. 



2 Finite element analysis of periodic structures 

Consider a periodic structure, as shown in figure 2, which is made of a large 
number N of cells. We are interested by the computation of the frequency 
response function for a point force excitation F = 1 somewhere in the struc- 
ture and a response u at another point. We propose here an efficient method 
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to compute this function by using a recursive approach to get the dynamic 
stiffness matrix of different sets of cells. 



2. 1 Behaviour of a cell 

Consider first the case of only one cell. The discrete dynamic equation of a 
cell obtained from a FE model at a frequency uj and for the time dependence 
g-iwt jg given by: 

(K - ia;C - a;2M)q = f (1) 

where K, M and C arc the stiffness, mass and damping matrices, respectively, 
f is the loading vector and q the vector of the degrees of freedom (dofs). A 
viscous damping is considered here but the same results could be obtained 
with other damping models. Introducing the dynamic stiffness matrix D = 
K — iLoC — a;^M, decomposing the dofs into boundary [B) and interior (7) 
dofs as shown in figure 3, and assuming that there are no external forces on 
the interior nodes, result in the following equation: 
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(2) 



The interior dofs can be eliminated using the second row of equation (2), which 
results in 

q/ = -Dy/D/sqs (3) 

The first row of equation (2) becomes 

fij = {pBB - Db/D7/D/b) qs (4) 

It should be noted that only boundary dofs are considered in the following. 
The cell is assumed to be meshed with an equal number of nodes on their 
opposite sides. The boundary dofs for one cell are decomposed into left (L) 
and right {R) dofs as shown in figure 3. Thus, equation (4) is rewritten as 
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where D^^^ is the dynamic stiffness matrix of a single ceU. This matrix is 
symmetric if the matrices K, M and C in relation (1) are symmetric. 



2.2 Computation of reduced dynamic stiffness matrices 



Consider a structure made of two cells with the respective dynamic stiffness 
matrices denoted by A and B as in figure 4. We propose to remove the internal 
degrees of freedom at the boundary between the two cells to compute the 

dynamic stiffness matrix, denoted D*^^\ relating the degrees of freedom (dofs) 
in the first section of A and the last section of B. The dynamic stiffness matrix 
of the substructure with two cells is computed by 



fl 




All 


Alr 







qi 
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Arr + Bll 
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qa 



As there is no load on the interior section, one gets f2 = and 

q2 = - {Arr + BllY^ {Arl^Ii + Bi^qs) (7) 



The global dynamic stiffness matrix of the two cells structure is thus 



fl 
fa 



All- Alr{Arr + Bll) ^ Arl -Alr{Arr + Bll) ^^lr 

— ^RL {Arr + B^l) ^ Arl Br/j — Brl {Arr + Bll) ^ B 



LR 





qi 




qa 



:D(2) 



qi 
qa 



This defines the matrix D^^^ which relates the forces and displacements dofs at 
the extreme sections of the two-cells structure. It can be easily checked that if 
the matrices A and B are symmetric, the resulting matrix D*^^) of relation (8) 
is also symmetric. The operation of removing the interior dofs is now denoted 
by {., .} such that we can write 



D(') = {A, B} 



(9) 
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2.3 Case of general structures 



Consider now a structure without internal load and made of 2" cells. We 
propose to recursively remove the internal degrees of freedom between adjacent 
cells to compute the dynamic stiffness matrix, denoted 0*^^"^ relating the 
degrees of freedom in sections 1 and 2". Consider firstly a structure with 
two identical cells. From the precedent analysis, one sees that its dynamic 
stiffness matrix is given by D^^^ = {D^^^ D^^^j. Repeating the process (see an 
illustration in figure 5) , one gets the dynamic stiffness matrix of the structure 
with 2" cells in n steps by the recursive relation 

D{2")^{d(2"-^),D(2"")} (10) 



This matrix is such that 



fl 




f2" 







qi 




qi 




q2" 




q2" 



In cases where the structure is not composed of a number of cells which equals 
a power of two, one can modify the previous procedure using the binary 
representation of the total number of cells N. Consider the example where 
= 11 = lOllfe in binary representation. One calculates first the dynamic 
stiffness matrix for a structure with 8 cells, then this structure is assembled 
with a structure made of 2 cells which has been computed during the compu- 
tation of the 8 cells structure. Finally the resulting matrix is assembled with 
a one cell matrix. The approach can be resumed by 

D(") = {{d(8),D(2)},D«} (12) 



in which the matrices D*^^") are computed by the method presented before. 
The elimination of the interior degrees of freedom gives the matrix linking 
the forces and displacements degrees of freedom in sections 1 and A^. One can 
notice that the final matrix is computed with a number of operations of order 
log2 N, thus saving a huge number of computations when we compare to a 
standard approach in which all the matrices of the cells are assembled into a 
global matrix. 

Using this method it is easy to compute the matrices for the parts of the 
structure respectively on the left and on the right of the force. Assembling 
these two matrices, applying the appropriate boundary conditions on the first 
and last sections, one gets the final linear system. This system has a number of 
dofs which equals approximately three times the number of dofs in a section. 
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3 Examples 



3.1 Beam structure 



In the first example, we consider the beam shown in figure 6 which is made 
of elements with four dofs. The stiffness and mass matrices of an element of 
length I are given by 



EI 



12 6Z -12 6Z 

-12 -6/ 12 -6/ 

6/ 2/2 -6/ 4/2 



(13) 



pSl 
420 



156 22/ 54 -13/ 
22/ 4/2 13/ -3/2 
54 13/ 156 -22/ 
-13/ -3/2 -22/ 4/2 



(14) 



Here, E is the Young's modulus, / the second moment of area, p the density 
of the material and S the cross-sectional area of the beam. Using the previous 
approach, one can compute the dynamic stiffness matrices for the sections on 
the left and on the right of the force. The damping matrix is obtained by using 
a complex Young modulus such that E = Eo^l + ir]) with r] = 0.01 leading to a 
hysteretic damping instead of the viscous one such that D = (l-|-i77)K — a;2M. 
The beam is made of steal such that Eq ^ 2 x lO^^Pa, p = 7800kg/m^, the 
length of the beam L = Im, / = 8.33 x 10~^^m^ and S = 10~^m2. The matrix 
for the whole beam is obtained by assembling the matrices of the left and right 
parts of the beam on each sides of the force. Taking into account the fixed 
displacement boundary conditions by removing the corresponding degrees of 
freedom in the global matrix results in the final system. Figure 7 presents 
the frequency response functions for structures with respectively 8 and 1024 
elements in each part of the beam. It can be seen that 8 elements are not 
sufficient to compute accurately the solution while 1024 elements lead to a 
very good result. 

In figure 8, the cpu time of the computation is plotted versus the global number 
of elements in the beam. The computation of 10000 points in frequency is made 
for each mesh of the beam and the largest mesh has 2 x 4196 elements for the 
complete beam. A linear behavior of the cpu time versus the logarithm of the 
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number of elements can be seen as expected. 



3.2 Plate 

Consider now the plate shown in figure 9. The mesh of a cell is obtained 
by Abaqus and consists in 50 elements of size Ly/50 x Lx/T^ where n is 
the number of cells along the direction x. The mass and stiffness matrices 
are produced by Abaqus then there are loaded in Matlab and the precedent 
procedure allows computing the displacement for a load at the centre of the 
plate. Information on the size are given in figure 9 and the plate is still made 
of steal. The boundary conditions are simply supported on all sides. Figure 10 
presents the frequency response functions for structures with globally 32 and 
512 cells. It can be seen, as for the beam, that 32 cells are not sufficient to 
compute accurately the solution while 512 cells lead to much better results. For 
high frequencies a discrepancy with the analytical solution can still be seen. 
It has been checked that the result can be considerably improved by taking 
100 elements instead of 50 along the direction y. The figures also present the 
results from the standard finite element approach obtained by assembling the 
elementary matrices for each period and solving the global linear system. It 
can be seen that both finite elements computations yield identical results. 

In figure 11, the cpu time of the computation is plotted versus the global 
number of cells in the plate. The computation of 100 points in frequency is 
made for each mesh of the plate and the largest mesh has 2048 cells along x 
for the complete plate. Once again a linear behavior of the cpu time versus 
the logarithm of the number of cells can be observed. The computing times 
for the recursive and standard finite element methods are compared in table 
1. It can be seen that the recursive method is a little slower than the standard 
method for a number of periods lower than 16. For a larger number of periods 
the recursive method tends to be more and more efficient as the number of 
periods increases. 



4 Conclusion 

A method has been described to compute frequency response functions for 
waveguide structures with periodic or homogeneous sections. The proposed 
method allows computing the solution in a time proportional to the logarithm 
of the numbers of cells in the structure. This simple recursive method can be 
applied for any type of one- dimensional periodic waveguides. It could be used 
in the future for the computation of complex structures such as tyres. 
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Fig. 1. Examples of waveguide structures. 
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Fig. 2. Periodic structure made of N cells with an excitation by a force F. 
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Fig. 3. Interior and boundary dofs for a single cell. 
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Fig. 4. Structure with two cells and three sections. 
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Fig. 5. Structure at each step. 
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Fig. 6. Beam structure (a) and beam element (b). 
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Fig. 7. Frequency response functions for a beam with 2x8 elements (upper graph) 

and 2 x 1024 elements (lower graph): analytical solution, finite element 

solution. The position of the excitation point is shown in figure 6 and the response 
is computed at the same point. 
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Number of elements 



;. 8. Computing time versus the number of cells in the beam. 



19 



L=lm 




X 



Fig. 9. Plate example. 
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Fig. 10. Frequency response functions for a plate with 2 x 16 cells (upper graph) and 

2 X 256 cells (lower graph): analytical solution, recursive finite element 

solution and . . standard finite element solution. The excitation is located at the 
center of the plate as shown in figure 9 and the response is computed at the same 
point. 
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Number of elements 



Fig. 11. Computing time versus the number of cells in the plate. 



22 



Table 1 



Computing times for recursive and standard FEM methods. 



number of periods 


recursive method (s) 


classical fem (s) 


2 


18.3 


14.4 


4 


28.8 


20.1 


8 


39.5 


33.7 


16 


50.3 


69.2 


32 


57.0 


164.3 


64 


71.8 


435.0 


128 


82.7 


1295.0 


256 


93.4 


4093.1 


512 


100.0 


13901.4 


1024 


114.8 


49271.1 


2048 


125.6 


184397.3 
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